Little NASA (Group 9): Vedant Srivastava, Weibo Chen, Xinyu Zhao¶
Introduction¶
Have you ever looked up at the sky and wondered, is our planet Earth the only place with life on it? The answer to that question is, probably not. With this arises the other question, which of the planets out there can actually support life on it? That time is not far when humans too might be able to / have to start living on planets other than Earth.
With this comes the significance of this project. Back in 2009, NASA launched The Kepler Space Observatory, a NASA-built satellite for space observation. As mentioned by NASA, the satellite is dedicated to finding exoplanets or extrasolar planet which by definition is a planet outside the Solar System. The ultimate goal of this exploration is the posibility of finding other habitable planets, besides our own. The telescope within the satellite has remained in service since 2014 on an extended mission and NASA consistently releases the data of it for public research interest.
We took this data as our working dataset of the project. Through manipulation of this data, we somewhat played the roles of amateur astronomers, in hopes of finding out which of the key components identify an astral body as an exoplanet and then predict the status of such other astral bodies based on what we have learnt.
Research Interests¶
Our project research questions are devided into 5 main parts:
- How is the proportion of exoplanets being confirmed in the existing literature but disconfirmed by measurements from Kepler? How about the other way round?
- How is the koi_score determined? Which parameters contribute the koi_score generation? Will the Transit Properties weight more in the determination of koi_score than the Stellar Parameters? Currently, we suppose the key predictors are koi_period, koi_impact, koi_duration, koi_prad, koi_teq, koi_smass and koi_sage, will that be the case?
- What regression models can we create? Which regression model will be the best chioce?
- Are these models really good for the explanatory and prediction?
- Still there are some KOI(Kepler Objects of Interest, astral bodies in Kepler mission) left undetermined, can we use our prediction model to make a reasonable guess of a given KOI being an exoplanet or not?
The objectives for our project will be:
- To see the proportion of the matched and mismatched classification of celestial bodies between previous archives and by Kepler Data Analysis score.
- To determine which parameters contribute the determination of a celestial body as an exoplanet or not based on literature review. In detail, we would like to find out how koi_score was determined by parametes from transit properties and stellar properties.
- To run the descriptive statistics for our data and summarize the characteristics of our variables.
- To use histogram to visualize numeric variables in our datasets.
- To build regression models.
- To visualize the model diagnostics.
- To split data into two parts, one with determined koi_score (model building set), one without (to-predict set). We will build and improve our prediction model from the model building dataset and then apply the model into the to-predict dataset to see in which category a certain celestial body falls.
Data¶
Section 1: Data Dictionary¶
Variables are separated into different sections for better understanding
Descriptions are captured from NASA Exoplanet Science Institute
https://exoplanetarchive.ipac.caltech.edu/docs/API_kepcandidate_columns.html
Our focus is to build a model for classifying an astral body. We decide to use the following three categories of variables to build our predicting model:
Identification Variable
- kepid: Target identification number, as listed in the Kepler Input Catalog (KIC).
- kepoi_name: A number used to identify and track a Kepler Object of Interest (KOI).
- kepler_name: Kepler number name in the form "Kepler-N," plus a lower-case letter, identifying the planet.
- koi_disposition: The category of this KOI from the Exoplanet Archive. Current values are CANDIDATE, FALSE POSITIVE, NOT DISPOSITIONED or CONFIRMED.
- koi_pdisposition: The pipeline flag that designates the most probable physical explanation of the KOI. Typical values are FALSE POSITIVE, NOT DISPOSITIONED, and CANDIDATE.
- koi_score: A value between 0 and 1 that indicates the confidence in the KOI disposition. For CANDIDATEs, a higher value indicates more confidence in its disposition, while for FALSE POSITIVEs, a higher value indicates less confidence in that disposition.
- koi_fpflag_nt: Like Flag A KOI whose light curve is not consistent with that of a transiting planet.
- koi_fpflag_ss: A KOI that is observed to have a significant secondary event, transit shape, or out-of-eclipse variability, which indicates that the transit-like event is most likely caused by an eclipsing binary.
- koi_fpflag_co: The source of the signal is from a nearby star.
- koi_fpflag_ec: The KOI shares the same period and epoch as another object and is judged to be the result of flux contamination in the aperture or electronic crosstalk.
- koi_disp_prov: Disposition Provenance.
Transit Properties
- koi_period: The interval between consecutive planetary transits.
- koi_time0bk: The time corresponding to the center of the first detected transit in Barycentric Julian Day (BJD) minus a constant offset of 2,454,833.0 days.
- koi_time0: The time corresponding to the center of the first detected transit in Barycentric Julian Day (BJD).
- koi_eccen: Eccentricity Value.
- koi_longp: Longitude of Periastron.
- koi_impact: The sky-projected distance between the center of the stellar disc and the center of the planet disc at conjunction.
- koi_duration:The duration of the observed transits.
- koi_ingress: The time between first and second contact of the planetary transit.
- koi_depth: The fraction of stellar flux lost at the minimum of the planetary transit.
- koi_ror: The planet radius divided by the stellar radius.
- koi_srho: Fitted stellar density is a direct observable from the light curve that, in the small-planet approximation, depends only on the transit's period, depth, and duration.
- koi_fittype: Type of Fit for planetary parameters.
- none
- LS: Least Squares fit
- MCMC: Markov Chain Monte Carlo fit
- LS+MCMC: Least Squares Fit with Markov Monte Carlo error bars
- DV: Data Validation pipeline fit
- koi_prad: The radius of the planet.
- koi_sma: Half of the long axis of the ellipse defining a planet's orbit.
- koi_incl: The angle between the plane of the sky (perpendicular to the line of sight) and the orbital plane of the planet candidate.
- koi_teq: Approximation for the temperature of the planet.
- koi_insol: Insolation flux is another way to give the equilibrium temperature.
- koi_dor: The distance between the planet and the star at mid-transit divided by the stellar radius.
Stellar Parameters
- koi_steff: The photospheric temperature of the star.
- koi_slogg: The base-10 logarithm of the acceleration due to gravity at the surface of the star.
- koi_smet: The base-10 logarithm of the Fe to H ratio at the surface of the star, normalized by the solar Fe to H ratio.
- koi_srad: The photospheric radius of the star.
- koi_smass: The mass of the star.
- koi_sage: The age of the star.
- koi_sparprov: A flag describing the source of the stellar parameters.
Determination of variables we interested¶
Based on our literature reivew, we finalized our variables as the followings:
Dependent Variable:
- koi_score: A value between 0 and 1 that indicates the confidence in the KOI disposition. For CANDIDATEs, a higher value indicates more confidence in its disposition, while for FALSE POSITIVEs, a higher value indicates less confidence in that disposition.
Independent Variables:
- koi_period: The interval between consecutive planetary transits.
- koi_duration:The duration of the observed transits.
- koi_depth: The fraction of stellar flux lost at the minimum of the planetary transit.
- koi_srho: Fitted stellar density is a direct observable from the light curve that, in the small-planet approximation, depends only on the transit's period, depth, and duration.
- koi_fittype: Type of Fit for planetary parameters.
- koi_prad: The radius of the planet.
- koi_sma: Half of the long axis of the ellipse defining a planet's orbit.
- koi_incl: The angle between the plane of the sky (perpendicular to the line of sight) and the orbital plane of the planet candidate.
- koi_teq: Approximation for the temperature of the planet.
- koi_dor: The distance between the planet and the star at mid-transit divided by the stellar radius.
- koi_steff: The photospheric temperature of the star.
- koi_slogg: The base-10 logarithm of the acceleration due to gravity at the surface of the star.
- koi_smet: The base-10 logarithm of the Fe to H ratio at the surface of the star, normalized by the solar Fe to H ratio.
- koi_srad: The photospheric radius of the star.
- koi_smass: The mass of the star.
Section 2: Data Cleaning¶
The data is downloaded from NASA Exoplanet Archive webpage:(https://exoplanetarchive.ipac.caltech.edu/docs/data.html). The data comes in a well-formatted, decent form, no heavy data cleaning work is required. We manipulate the dataset based on our own research interests: filtering out irrelavent variables, dropping the NULL data input, assigning the correct column type for each column, creating dummy variables for categorical variables (dummy variables creating is done in regression part).
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
nasarawdata=pd.read_csv('Nasa.csv')
nasarawdata.shape
Variable Selection¶
regvarlist=[#Project Disposition Columns
'koi_score',
#Transit, and Planet Properties
'koi_period','koi_duration','koi_depth',
'koi_srho','koi_fittype','koi_prad','koi_sma','koi_incl','koi_teq','koi_dor',
#Stellar Parameters
'koi_steff','koi_slogg','koi_smet','koi_srad','koi_smass']
nasadataregvar=nasarawdata[regvarlist].copy()
nasadataregvar.dtypes
All columns are read in correct types.
#Delete observations contain any null values
nasadatareg=nasadataregvar.dropna()
nasadatareg.shape
In total, we have 7989 vaild oberstavions. We have 1 dependent variable and 15 independent variables.
Data Analysis¶
Objective 1: Literature & KOI concordance analysis¶
From prelinimary literature review, we got the meaning of all the columns so as to solve the first the question: How is the proportion of exoplanets being confirmed in the existing literature but disconfirmed by measurements from Kepler? How about the other way round?
We used raw data to conduct this analysis.
kepmatchstat=nasarawdata[['kepid','koi_disposition','koi_pdisposition']].copy()
kepmatchstat.groupby(["koi_disposition", "koi_pdisposition"]).size()
Because in Kepler mission, to be conservative, no 'CONFIRMED' value for koi_pdisposition column, we would treat observations having 'CONFRIMED' in koi_disposition and 'FALSE POSITIVE' in koi_pdisposition, and observations having 'FALSE POSITIVE' in koi_disposition and 'CANDIDATE' in koi_pdisposition as mismatching.
The python output indicates that there are in total 8 mismatching cases out of 9564 observations, an extreme similarity between archive-determined exoplanets and Kepler exoplanet determining algorithm. The result brings us confidence to use literature derived key parameters to construct (or deduct) the algorithm that the Kepler mission might use.
Objective 2 Predicting Model¶
Section 1 Linear regression modeling¶
2.1.1 Descriptive Statistics¶
nasadatareg.head(5)
Data description¶
nasadatareg.describe().T
The descriptive statstics give us a brief understanding of our data. The percentile data tells us that our data have heavy skewness and extreme values. We will explore more insights in the following analyses.
Independent variable visualizations¶
hisvarlist=['koi_period','koi_duration','koi_depth','koi_srho','koi_prad','koi_sma','koi_incl','koi_teq','koi_dor',
'koi_steff','koi_slogg','koi_smet','koi_srad','koi_smass']
for i in hisvarlist:
plt.figure(figsize=(8,5))
plt.title("{}".format(i),fontsize=16)
plt.hist(nasadatareg[i],color='cyan',edgecolor='k')
plt.grid(True)
plt.xlabel(i,fontsize=14)
plt.show()
We visualized each column (only continuous one) by plotting the histogram, we found that most of the variables are heavily skewed.
- koi_period, koi_duration, koi_depth, koi_srho, koi_fittype, koi_prad, koi_sma, koi_teq, koi_dor, koi_steff, koi_srad, koi_smass are right-skewed
- koi_incl, koi_slogg, koi_smet are left-skewed
Such visualizations gave us an impression that for our modeling purpose, perhaps applying certain transformations on the indepedent variables is better than using the identical form.
Dependent variable visualization¶
plt.figure(figsize=(8,5))
plt.title("koi_score",fontsize=16)
plt.hist(nasadatareg['koi_score'],color='cyan',edgecolor='k')
plt.grid(True)
plt.xlabel("koi_score",fontsize=14)
plt.show()
The vast majority of the koi_score is 0 and 1.
The violation of normal distribution indicates that a linear regression model might not be that robust.
Moreover, the nearly dichotomized koi_score indicates a logistic regression model could be a good alternative to serve our purpose.
Visualizing the association between Dependent variable and Independent variables¶
for c in nasadatareg.columns[1:]:
plt.figure(figsize=(8,5))
plt.title("{} vs. \nkoi_score".format(c),fontsize=16)
plt.scatter(x=nasadatareg[c],y=nasadatareg['koi_score'],color='blue',edgecolor='k')
plt.grid(True)
plt.xlabel(c,fontsize=14)
plt.ylabel('koi_score',fontsize=14)
plt.show()
For transit property plots (plots before _koi_steff vs. koiscore ), most of them are in 'C' shape, showing plenty of observations clustered around the upper limit (koi_score=1) and the lower limit (koi_score=0).
For stellar parameter plots (plots after _koi_steff vs. koiscore ), most of them are in 'I' shape. Same as transit property plots, lots of the observations are clustered around the upper limit and the lower limit; the difference is the observations are more scattered around than in transit property plots.
Also, the scatterplots show some outliers,
- koi_period vs. koi_score, bottom right corner
- koi_duration vs. koi_score, bottom right corner
- koi_srho vs. koi_score, bottom right corner and upper right corner
- koi_prad vs. koi_score, bottom right corner
- koi_teq vs. koi_score, bottom right corner and upper right corner
- koi_dor vs. koi_score, bottom right corner
- koi_steff vs. koi_score, bottom right corner
- koi_slogg vs. koi_score, upper left corner
- koi_srad vs. koi_score, bottom right corner
These are not good signs for linear regression because from the scatterplots hardly can we find a linear trend line to describe the relationship of koi_score along with each paramter (a perpendicular line to the X-Axis is not a feasible trend), and our linear model may bias becasue of the outliers. However, the clustering feature is a good sign for applying logistic regression.
Create dummy variables for category variable¶
nasadatareg2=pd.get_dummies(nasadatareg)
nasadatareg3=nasadatareg2.rename(columns={"koi_fittype_LS+MCMC": "koi_fittype_LS_MCMC"})
Correlation matrix¶
We also want to know about the associations between Dependent variables
nasadatareg3.head()
corr = nasadatareg3[nasadatareg3.columns[:-2]].corr()
corr
Correlation matrix is a basic method to view pairwise correlation coefficients. However, it is a little bit hard for distinguishing which pairs are highly correlated. We decide to use the pairwise scatter plots and heatmap for more informative displays.
Variable pairwised scatter plots¶
from seaborn import pairplot
pairplot(nasadatareg3[nasadatareg3.columns[:-2]])
We generate a visualized correlation matrix in scatter plot using pairplot in Seaborn. From this plot matrix we can find some correlated pairs:
- koi_sma vs. koi_period
- koi_srad vs. koi_slogg
- koi_smass vs. koi_slogg
- koi_steff vs. koi_slogg
Heatmap method¶
from statsmodels.graphics.correlation import plot_corr
sns.set(rc={'figure.figsize':(12,7)})
fig=plot_corr(corr,xnames=corr.columns,title='Correlation Heatmap')
The heatmap outstandingly helps us to find out which pairs are highly correlated. The girds have deep colors (Both in deep red and in deep blue) indicates the corresponding pairs are highly correlated.
The highly correlated pairs are:
- koi_sma vs. koi_period
- koi_dor vs. koi_period
- koi_dor vs. koi_sma
- koi_slogg vs. koi_srad
- koi_slogg vs. koi_smass
We decide to keep all of them into the saturated model although we know that will introduce multicollinearity issue. The reason is that we consider those parameters are crucial based on our literature review.
2.1.2 Modeling¶
Library statsmodels is used to generate and test our predicting model.
import statsmodels.api as sm
import statsmodels.formula.api as smf
At the beginning, we need to specify the linear model, a '~' sign represents the equal sign connecting the dependant variable and independant variables.
model_str = nasadatareg3.columns[0]+' ~ '+'+'.join(nasadatareg3.columns[1:-2])+'+'+nasadatareg3.columns[-1]
model_str
model=smf.ols(formula=model_str, data=nasadatareg3)
fitted = model.fit()
print(fitted.summary())
We select koi_teq and koi_fittype_LS_MCMC as examples to interpret our model.
- koi_teq
On average, for each Kelvin unit increase in planet temperature while holding other variables constant, the koi_score will be decreased by 0.0001
- koi_fittype_LS_MCMC
Hypothetically, on average, for an identical astral body, LS+MCMC fitting method will have 0.1678 higher score than LS fitting method.
result=pd.DataFrame()
result['pvalues']=fitted.pvalues[1:]
paramlist=nasadatareg3.columns[1:-2].tolist()
paramlist.append('koi_fittype')
result['Parameters']=paramlist
result.set_index('Parameters',inplace=True)
result['Statistically significant?']= result['pvalues'].apply(lambda x:'Yes' if x<0.05 else 'No')
result
R-square is 0.341 and Adj. R-squared is 0.34, both mean our predicting model is not that capable of doing the predition work since more than half dependant variance is left unexplained.
Parameters koi_srho, koi_prad, koi_dor, koi_steff are not statistically significant at 0.05 level, however, since we valued literature reivewed model more, and p-values are acceptable, we decide to keep them in the model.
Section 2 Model Diagnosis¶
We have our model now, but we want to know more about how well, and robust our model is. This sector is mainly about residual analysis to check if certain linear regression assumptions are held or not.
2.2.1 Residual Pattern: Homoscedasticity OR Heteroscedasticity?¶
Residuals vs. predicting variables plots¶
for d in nasadatareg[hisvarlist]:
plt.figure(figsize=(8,5))
plt.title("{} vs. \nModel residuals".format(d),fontsize=16)
plt.scatter(x=nasadatareg[d],y=fitted.resid,color='blue',edgecolor='k')
plt.grid(True)
xmin=min(nasadatareg[d])
xmax = max(nasadatareg[d])
plt.hlines(y=0,xmin=xmin*0.9,xmax=xmax*1.1,color='red',linestyle='--',lw=3)
plt.xlabel(d,fontsize=14)
plt.ylabel('Residuals',fontsize=14)
plt.show()
Residual plots show some bits of clustering. Some of them have an oblique line at the bottom, suggesting the residuals follow a slightly positive trend. These characteristics tell us that the residual are following some undesired pattern, and thus the assumption of homoscedasticity is violated.
Also, this non-random pattern indicates that the deterministic portion of this model is not capturing the whole explanatory information that is affecting the residuals. Possibilities that could cause such patterns include:
- A missing variable
- A missing higher-order term of a variable in the model to explain the curvature
- A missing interaction between terms already in the model
Some of the plots also show residual outliers that are away from the main cluster.
Fitted vs. residuals¶
plt.figure(figsize=(8,5))
p=plt.scatter(x=fitted.fittedvalues,y=fitted.resid,edgecolor='k')
xmin=min(fitted.fittedvalues)
xmax = max(fitted.fittedvalues)
plt.hlines(y=0,xmin=xmin*0.9,xmax=xmax*1.1,color='red',linestyle='--',lw=3)
plt.xlabel("Fitted values",fontsize=15)
plt.ylabel("Residuals",fontsize=15)
plt.title("Fitted vs. residuals plot",fontsize=18)
plt.grid(True)
plt.show()
This plot looks pretty bad. Again, a clear negative trend of residuals vs. fitted values indicates that we must have omitted one or more variables from our model that could be account for this trend line. However, we couldn't know exactly which one(s) we missed.
2.2.2 Normality of residuals¶
We applied both raw residuals and standardized residuals for normality test.
1. Histogram of standardized residuals¶
plt.figure(figsize=(8,5))
plt.hist(fitted.resid_pearson,bins=20,edgecolor='k')
plt.ylabel('Count',fontsize=15)
plt.xlabel('Standardized residuals',fontsize=15)
plt.title("Histogram of standardized residuals",fontsize=18)
plt.show()
The histogram shows that the normalized residuals are neither symmetric nor bell-shape distributed.
2. Q-Q plot of the residuals¶
from statsmodels.graphics.gofplots import qqplot
plt.figure(figsize=(8,5))
fig=qqplot(fitted.resid_pearson,line='45',fit='True')
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
plt.xlabel("Theoretical quantiles",fontsize=15)
plt.ylabel("Sample quantiles",fontsize=15)
plt.title("Q-Q plot of standardized residuals",fontsize=18)
plt.grid(True)
plt.show()
The Q-Q plot shows that, although not perfectly aligned, the standardized residuals are generally following the diagnol reference line.
3. Normality (Shapiro-Wilk) test of the residuals¶
from scipy.stats import shapiro
_,p=shapiro(fitted.resid)
if p<0.01:
print("The residuals seem to follow normal distribution")
else:
print("The normality assumption may not hold")
Finally we used Shapiro-Wilk test for testing the normality. The test result shows the residuals are very likely following normal distribution.
2.2.3 Outliers detection¶
Cook's distance¶
from statsmodels.stats.outliers_influence import OLSInfluence as influence
inf=influence(fitted)
(c, p) = inf.cooks_distance
plt.figure(figsize=(8,5))
plt.title("Cook's distance plot for the residuals",fontsize=16)
plt.stem(np.arange(len(c)), c, markerfmt=",")
plt.grid(True)
plt.show()
Cook's distance is for outliers detection. As the plot shows, there is one significant outlier in our dataframe. In our future model improving process we will remove this outlier and rerun the model to see how well our model could be improved.
2.2.4 Multicollinearity Detection¶
Variance inflation factor¶
from statsmodels.stats.outliers_influence import variance_inflation_factor as vif
for i in range(len(nasadatareg3.columns[1:])):
v=vif(np.matrix(nasadatareg3[1:]),i)
print("Variance inflation factor for {}: {}".format(nasadatareg3.columns[i],round(v,2)))
There are two parameters with VIF > 10, thereby indicating significant multicollinearity. Because we intentionally included highly correlated exogenous variables in our model, this result is not beyond our expection.
2.2.5 Linear Model & Diagnosis Conclusion¶
We successfully built the model for predicting the koi_score of an astral body according to its own properties. However this model has several flaws that could be corrected in our future analysis.
- The heavily skewed distribution of each independent variables suggestes proper transformation may improve our model. This is also logically reasonable since the physcial laws applied for deriving the parameters of astral body require heavily transformations.
- Apart from the parameters we chose, we probably missed some key parameters in our model becasue
- The Adj R-squared is quite low.
- Most of the residual plots follow some patterns.
- There are outliers in our model, our model might be more stable (p-value decrease) after the outliers are removed from our dataset.
Section 3: Logistic Regression modeling under Statsmodel¶
Recall the pattern of the histogram of dependent variable (koi_score), vast majority of the observations are 0 and 1. We consider Logistic Regression might perform better than linear regression model.
2.3.1 Descriptive Statistics¶
nasadatalogvar=nasarawdata[regvarlist].copy()
For logistic regression, the outcomes are binary values (e.g. 1 or 0, Yes or No). Therefore, we dichotomize the koi_score into 0 or 1 according the criteria: koi_score value above 0.5 will be assigned value 1, otherwise it will be assigned value 0.
#Dichotomize koi_score
nasadatalogvar['score_bi']=nasadatalogvar['koi_score'].apply(lambda x:0 if x<=0.5 else 1)
#Row-wise drop null values
nasadatalog=nasadatalogvar.dropna()
nasadatalog['score_bi'].value_counts()
sns.countplot(x='score_bi', data=nasadatalog, palette='hls')
plt.show()
In total, we have 7989 observations with complete observation values.
The two groups looks well balanced, there is no need for oversampling to balance the groups.
meancomparison=nasadatalog.groupby('score_bi').mean()
for i in meancomparison.columns[1:]:
plt.figure(figsize=(8,5))
plt.title("{} vs. \nscore_bi".format(i),fontsize=16)
plt.bar(meancomparison.index,meancomparison[i],color='blue',edgecolor='k')
plt.grid(True)
plt.ylabel(i,fontsize=14)
plt.xlabel('score_bi',fontsize=14)
plt.show()
We plotted the bar charts of two score group means to do a naked eye evalutation of, for each variable, the existance of distinguishable difference between two groups.
Considering the magnitutde of Y-axis, we can only conclude that the mean value of koi_depth in two score groups are substantially different.
For more statistical conclusions, we will rely on two sample t-test to determine if there exists statistical difference between two groups for each independent variable.
Two sample t-test¶
Before applying the two sample t-test, we have to determine if the variances of the two groups are equal.
Recall the histogram of each dependent variable, most of them are heavily skewed, Levene Test is robust for this situation.
The process is
- Applying Levene Test to determine the equality of variance.
- Determine the equal_var argument based on the Levene Test result and perform preferred tests.
import scipy.stats as ss
from scipy.stats import levene
nasadatalog.groupby('score_bi').mean()
twos_ttest=[]
for i in hisvarlist:
#Variables 'koi_steff','koi_smet' and 'koi_smass' are not so heavily skewed, set center argument to 'median'
if i in ['koi_steff','koi_smet','koi_smass']:
if levene(nasadatalog[nasadatalog['score_bi']==1][i],
nasadatalog[nasadatalog['score_bi']==0][i],center='median').pvalue<0.05:
#If p-value is less than 0.05, variances of two groups are considered unequal
twos_ttest.append([i,ss.ttest_ind(a=nasadatalog[nasadatalog['score_bi']==1][i],
b=nasadatalog[nasadatalog['score_bi']==0][i], equal_var=False).pvalue])
else:
#Otherwise two groups are considered equal
twos_ttest.append([i,ss.ttest_ind(a=nasadatalog[nasadatalog['score_bi']==1][i],
b=nasadatalog[nasadatalog['score_bi']==0][i], equal_var=True).pvalue])
#Except those three variables, the rest variables are heavily skewed, set center argument to 'trimmed'
else:
if levene(nasadatalog[nasadatalog['score_bi']==1][i],
nasadatalog[nasadatalog['score_bi']==0][i],center='trimmed').pvalue<0.05:
twos_ttest.append([i,ss.ttest_ind(a=nasadatalog[nasadatalog['score_bi']==1][i],
b=nasadatalog[nasadatalog['score_bi']==0][i], equal_var=False).pvalue])
else:
twos_ttest.append([i,ss.ttest_ind(a=nasadatalog[nasadatalog['score_bi']==1][i],
b=nasadatalog[nasadatalog['score_bi']==0][i], equal_var=True).pvalue])
df_twos_ttest=pd.DataFrame(twos_ttest,columns=['Parameter','p-value'])
df_twos_ttest['p<0.05']=df_twos_ttest['p-value'].apply(lambda x:'Yes' if x <0.05 else 'No')
df_twos_ttest.set_index('Parameter')
koi_srho, the stellar density, is the only one that no statistically significant difference exists in two score groups.
Similar to what we've done in linear model, we value the literature reivewed model more and this variable will still be included into our saturated model.
2.3.2 Modeling¶
We use logit model under statsmodel, discrete_model library for modeling logistic regression.
import statsmodels.discrete.discrete_model as smdd
nasadatalog2=pd.get_dummies(nasadatalog)
varlist=['koi_period','koi_duration','koi_depth','koi_srho','koi_prad','koi_sma','koi_incl','koi_teq','koi_dor',
'koi_steff','koi_slogg','koi_smet','koi_srad','koi_smass','koi_fittype_LS+MCMC']
X=sm.add_constant(nasadatalog2[varlist])
y=nasadatalog2['score_bi']
logit_model1=smdd.Logit(y,X)
result=logit_model1.fit()
print(result.summary2())
To our surprise, the p-value for koi_srho is 0.9306, which is too great and meaning this variable makes almost no contribution for the outcome classification. We decide to remove it from the saturated model and rerun the reduced model.
varlist2=['koi_period','koi_duration','koi_depth','koi_prad','koi_sma','koi_incl','koi_teq','koi_dor',
'koi_steff','koi_slogg','koi_smet','koi_srad','koi_smass','koi_fittype_LS+MCMC']
X2=sm.add_constant(nasadatalog2[varlist2])
y2=nasadatalog2['score_bi']
logit_model2=smdd.Logit(y2,X2)
result=logit_model2.fit()
print(result.summary2())
The AIC decreased and Log-Likelihood didn't change, meaning this reduced model is better than the saturated model. Although koi_sma is not significant, we consider this p-value level is accpetable and keep it as our finalized our model.
Using koi_steff and koi_fittype as two examples,the interpretations of our model are as follow.
- koi_teq
On average, for each Kelvin unit increase in planet temperature while holding other variables constant, the log odds ratio of koi_score=1 over koi_score=0 will be decreased by 0.0018.
- koi_fittype_LS+MCMC
Hypothetically, on average, for an identical astral body, LS+MCMC fitting method will have 1.6369 unit higher log odds ratio of koi_score=1 over koi_score=0 than LS fitting method.
2.3.3 Logistic Model Conclusion¶
- The logistic regression model has some improvements when compared to the linear regression model:
- Imporved Adj R-square.
- Decreased p-values.
- There is one major limitation in this model. The way we dichotomize the koi_score eliminates the magnitute of tendency. For example, A koi object having score 0.99 is different to a koi object having score 0.51, but due to our dichotomization, their intrinsic difference is eliminated.
- Similar to the linear regression model, the logistic model could be improved by:
- Applying transformation on independent variables.
- Adding the missing variables.
- Excluding the outliers.
Section 4: Machine Learning under SciKit Learn¶
We also want to know how well our logistic model is. One of the verification methods is to use machine learning technology and find out the prediction accruancy.
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
2.4.1 Prediction Accuracy¶
X_train, X_test, y_train, y_test = train_test_split(X2, y2, test_size=0.3, random_state=111111)
logreg = LogisticRegression(solver='liblinear')
logreg.fit(X_train, y_train)
y_pred = logreg.predict(X_test)
print('Accuracy of logistic regression classifier on test set: {:.2f}'.format(logreg.score(X_test, y_test)))
The predicting accuracy is way better than we thought, about 0.81.
2.4.2 ROC Curve¶
Another way to visualize how well the model in prediction is to plot the ROC Curve.
The receiver operating characteristic (ROC) curve is another common tool used with binary classifiers. The dotted line represents the ROC curve of a purely random classifier; a good classifier stays as far away from that line as possible (toward the top-left corner).
logit_roc_auc = roc_auc_score(y_test, logreg.predict(X_test))
fpr, tpr, thresholds = roc_curve(y_test, logreg.predict_proba(X_test)[:,1])
plt.figure()
plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()
Objective 3: Predicting the undetermined koi_score¶
The last objective is to use our regression model to predict the undetermined koi_score, through which, we can classify the kepler is a candidate or not.
#Finding the rows which have missing value for koi_score:
miss_score=nasarawdata[nasarawdata.koi_score.isnull()]
to_predict=miss_score[regvarlist].copy()
#Filtering out the rows that have missing value only for the koi_score column.
to_predict=to_predict[to_predict.isnull().sum(axis=1)<2]
to_predict.head()
We find something anormal, the koi_fittype is 'MCMC' type, which does not exist in our modeling dataset.
to_predict.koi_fittype.value_counts()
The worst case happen. The 'for prediction' dataset has one completely different feature to the modeling dataset: the koi_fittype that exists is only 'MCMC', and our model does not have a parameter coefficient for koi_fittype=MCMC. Therefore, we cannot use our model to predict the missing score because of the data structure mismatch.
Final Conclusions & Thoughts¶
In conclusion, we were surprised and satisfied with the well performance in prediction of our logistic model. For a given astral body with sufficient parameters, we could use this logistic model to predict in which category this astral body will likely to fall with 81% accuracy.
However, the logistic model itself, as well as the linear model, is still not good and requires improvements.
As we mentioned before, there are at least three major flaws wihtin our dataset and modeling procedure.
- Variable transformations
- Missing explanatory variable(s)
Outliers
We might generate more robust models after we addressed those issues.
Although we did heavy works on literature review, we found it is not enough for us to improve our model.
- We narrowed down the variable scope to the categories 'Transit Properties' and 'Stellar Parameters'. Depsite unselected variables within our pre-determined variable scope, the omitted variables may be outside this scope.
- Other categories describing other aspects of an astral body.
- Threshold-Crossing Event (TCE) Information
- KIC Parameters
- Pixel-Based KOI Vetting Statistics
- Other categories describing other aspects of an astral body.
- Understanding the terminologies of Cosmology and Astrophysics is challenging. The authors assume readers have knowledges about them, unfortunately we are not their target readers.
Figuring out proper transformations for our independent variables is associated with understanding the physical laws, and we don't have a member who is familiar with that.
To further improve our model, maybe we should ask help from a professional astronomer.
- We narrowed down the variable scope to the categories 'Transit Properties' and 'Stellar Parameters'. Depsite unselected variables within our pre-determined variable scope, the omitted variables may be outside this scope.
- We tried linear model (link function=identity) and logit model (link function=logit). We found logit link function is better than identity link function. However, there maybe other link functions that can perform much better than those two functions. In our future step, we can do box-cox transformation analysis (data-driven analysis) and take physical laws into consideration (theory-driven analysis) to find out a best-fit link function to connect the independent variables and dependent variable.
Reference¶
- Data source & Data dictionary
- NASA Exoplanet Archive
https://exoplanetarchive.ipac.caltech.edu/docs/data.html - Data Columns in Kepler Objects of Interest Table
https://exoplanetarchive.ipac.caltech.edu/docs/API_kepcandidate_columns.html
- NASA Exoplanet Archive
- Literature
- The NASA Exoplanet Archive: Data and Tools for Exoplanet Research
R. L. Akeson, X. Chen, D. Ciardi, M. Crane, J. Good, M. Harbut, E. Jackson, S. R. Kane, A. C. Laity, S. Leifer
https://iopscience.iop.org/article/10.1086/672273/meta - Exploring exoplanet populations with NASA’s Kepler Mission
Natalie M. Batalha
https://www.pnas.org/content/111/35/12647.full - Improved Parameters for Extrasolar Transiting Planets
Guillermo Torres, Joshua N. Winn, and Matthew J. Holman
https://iopscience.iop.org/article/10.1086/529429/meta - Improving Stellar and Planetary Parameters of Transiting Planet Systems: The Case of TrES-2
Alessandro Sozzetti, Guillermo Torres, David Charbonneau, David W. Latham, Matthew J. Holman, Joshua N. Winn, John B. Laird, and Francis T. O'Donovan
https://iopscience.iop.org/article/10.1086/519214/meta - The Transit Light Curve (TLC) Project. VI. Three Transits of the Exoplanet TrES-2
Matthew J. Holman, Joshua N. Winn, David W. Latham, Francis T. O'Donovan, David Charbonneau, Guillermo Torres, Alessandro Sozzetti, Jose Fernandez, and Mark E. Everett
https://iopscience.iop.org/article/10.1086/519077/meta
- The NASA Exoplanet Archive: Data and Tools for Exoplanet Research
- Python coding tutorials
- Machine-Learning-with-Python
tirthajyoti
https://github.com/tirthajyoti/Machine-Learning-with-Python/blob/master/Regression/Regression_Diagnostics.ipynb - Regression diagnostics
http://www.statsmodels.org/devel/examples/notebooks/generated/regression_diagnostics.html - Linear Regression Diagnostic in Python with StatsModels
Zhiya Zuo
https://zhiyzuo.github.io/Linear-Regression-Diagnostic-in-Python/ - Building A Logistic Regression in Python, Step by Step
Susan Li
https://towardsdatascience.com/building-a-logistic-regression-in-python-step-by-step-becd4d56c9c8 - Logistic Regression in Sci-Kit Learn
Stephen Barter
https://blog.goodaudience.com/classifying-flowers-using-logistic-regression-in-sci-kit-learn-38262416e4c6
- Machine-Learning-with-Python